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ABSTRACT 

Within a sufficiently large cosmic volume, conservation of baryons implies a simple ‘closed 
box’ view in which the sum of the baryonic components must equal a constant fraction of 
the total enclosed mass. We present evidence from Rhapsody-G hydrodynamic simulations 
of massive galaxy clusters that the closed-box expectation may hold to a surprising degree 
within the interior, non-linear regions of haloes. At a hxed halo mass, we hnd a significant 
anti-correlation between hot gas mass fraction and galaxy mass fraction (cold gas -I- stars), 
with a rank correlation coefficient of —0.69 within i? 5 ooc- Because of this anti-correlation, 
the total baryon mass serves as a low-scatter proxy for total cluster mass. The fractional 
scatter of total baryon fraction scales approximately as 0.02(Ac/100)°'®, while the scatter 
of either gas mass or stellar mass is larger in magnitude and declines more slowly with 
increasing radius. We discuss potential observational tests using cluster samples selected 
by optical and hot gas properties; the simulations suggest that joint selection on stellar 
and hot gas has potential to achieve 5% scatter in total halo mass. 

Key words: methods: numerical — galaxies: clusters: general — cosmology: theory — 
X-rays: galaxies: clusters. 


1 INTRODUCTION 


The abundance of galaxy clusters as a function of cluster mass 
is sensitive to both the growth of structure and cosmic expan¬ 
sion, providing not only stringent constraints on cosmological 
parameters but also consistency checks for the theory of gravity 
(see e.g. Miller et al. 2001[ Vikhlinin et al.|2009[ | Mantz et al. 


|2010||Rozo et al.||2010||Rapetti et al.||2013||Benson et al.||2013| 


and Allen et al.||2011 for a review). In cluster cosmology, the 

key is to accurately infer the mass of galaxy clusters from their 
observable properties, including gas mass and temperature 
from X-ray emission (e.g. Mantz et al.|2014 l, galaxy content 
from imaging and galaxy dynamics from spectroscopy (e.g. 
Kravtsov et al. 2014 Mamon et al. [20131, and strong and 


weak gravitational lensing effects (e.g. von der Linden et al. 


20141. Each of these mass proxies exhibits a certain amount of 


scatter around the true mass; minimizing and characterizing 
this scatter is essential for precision cosmology from galaxy 
cluster surveys (e.g. Lima &: Hu|2005 Wu et al.|2010 1. 


* Present address: California Institute of Technology, MC 367-17, 
Pasadena, CA 91125, USA. E-mail: hywu@caltech.edu 


To achieve accurate mass measurements, multi-wavelength 
observations have often been conducted for the same sample of 
clusters; for example, the CLASH project includes comparison 
between mass proxies from weak lensing. X-ray, and velocity 


dispersion | 

Postman et al.|2012 

Donahue et al.|2014 Biviano 

et al. 2013 

1; clusters observed by the South Pole Telescope 


using the Sunyaev-Zeldovich effect have been followed up pho¬ 
tometrically and spectroscopically ( Song et al.|2012 Ruel et al. 
20141. When multiple mass tracers are available for the same 


sample of galaxy clusters, a joint selection can reduce the mass 
scatter. In particular, the reduction of mass scatter is most 
effective when two mass tracers are anti-correlated with each 
other at a given mass (e.g. jCunhaj 2009| jStanek et al.||^10| 
Rozo et al.|2014 Evrard ^ al. 120141. 


Hydrodynamical simulations of galaxy clusters have been 


a powerful tool for understanding mass proxies (e.g. 

Evrard 

et al.j 1996 

Kravtsov et al.||2006| Rasia et al.| 2006 

Nagai 

et al.|2007| 

Stanek et al.|2010 Fabjan et al.|20111 [Rasia et al. 

2012 [Angulo et al.|2012[ Saro et al.||201 

B1 and the evolution 

of gas and stellar mass in clusters (e.g. 1 

Kravtsov et al.|2005 

Ettori et al.||2006| |Puchwein et al.||2010| |Young et al.||2011 
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Battaglia et al.|2013 Planelles et al.|2013 1. Recent results have 
shown that it is necessary to include the feedback of active 
galactic nuclei (AGN) in order to prevent catastrophic over¬ 
cooling in the cluster core, thus bringing the star formation in 
massive galaxies down to realistic levels and producing overall 
stellar mass fractions in better agreement with observations 


(e.g. Springel||2005 Sijacki et al.|2007 Booth fc Schaye|2009 


McCarthy et al.|20io Teyssier et al.|20TT Le Brun et al.|20T4 K 

In this work, we study gas and stellar mass fractions in 
a new set of hydrodynamical simulations of massive haloes. 
We find a significant anti-correlation between gas and stellar 
mass fractions that persists into the deeply non-linear regime. 
This anti-correlation does not simply reflect the well-known 
trends in the mean component fractions with total mass; the 
anti-correlation exists for deviations about the mean trends, 
meaning it reflects statistical behaviour of the component 
fractions at fixed halo mass. 

The new set of simulations is selected from the A^-body 
simulation sample Rhapsody (Wu et al. 2013a|b I and re¬ 
simulated with gas; we thus name our new sample Rhapsody- 
G. The original Rhapsody sample has been developed with 
the aim of understanding the impact of formation history on 
various mass tracers of galaxy clusters. In this paper, we focus 
only on the gas and stellar mass of the Rhapsody-G clusters; 
in companion papers (Hahn et al., in preparation and Martizzi 
et ah, in preparation), we will present detailed comparison 
between our simulations and observational results, including 
the properties of the BCG and the stellar mass-halo mass 
relation. 

This paper is organized as follows. Section introduces 
the simulations. In Section]^ we discuss the anti-correlation 
between gas and stellar mass fractions, while in Section]^ we 
discuss using total baryon mass as a low-scatter cluster mass 
proxy. We discuss the observational implications in Section 
and summarize our results in Section]^ 

Throughout this paper, we use radial and mass scales 
defined by a spherical density contrast with respect to the 
critical density of the Universe; e.g. Rsooc indicates the radius 
within which the average density is 500 times the critical density 
at the redshift of interest. 


2 SIMULATIONS 


The current Rhapsody-G simulation suite includes 10 hy¬ 
drodynamical zoom-in simulations centred on massive haloes 
from the original Rhapsody sample. Nine are chosen to have 
similar final mass of M 200 ~ 6 x and the tenth has 

M 200 ~ 1.3 X IO'^^Mq. We use the adaptive mesh refinement 


code Ramses (Teyssier 2002) and incorporate cooling, star 


formation, and AGN feedback. We describe here the simula¬ 
tion methods and recipes, as well as the details of our sample 
selection. 


2.1 Simulation methodology 


First, we briefly summarize the methods used to generate and 
post-process the simulations. We kindly refer the reader to 
Hahn et al. (in preparation) for more details. 

Precursor simulations. The galaxy clusters of the 
Rhapsody-G simulations are based on the Rhapsody A-body 


simulations ( Wu et al.|2013a b I, which include 96 cluster-sized 
haloes of mass Mvir = re-simulated with a 


mass resolution of 1.3 x lO®h“^M 0 . These haloes have been 
identified at z = 0 in a cosmological volume of 1 h~^Gpc^ 
from the LasDamas simulation suite. The 10 Rhapsody-G 
simulations presented here are selected from the full Rhap¬ 
sody sample in such a way that three of the main haloes have 
extreme concentration, two have an extreme number of sub¬ 
haloes, and five have approximately the median concentration 
and typical number of subhaloes. 


Initial conditions. We use the Music code I Hahn & Abel 


20111 to generate the initial conditions of the hydrodynamical 


simulations, at a starting redshift of 50. The Music code im¬ 
plements the second-order Lagrangian perturbation theory to 
generate displacements and velocities of dark matter particles, 
and the local Lagrangian approximation to generate a consis¬ 
tent initial density field of baryons on the grid. The Lagrangian 
volumes for the zoom simulations have been chosen to include 
a sphere of 8 /i“^Mpc centred on the clusters at 2 = 0, which 
allows us to study a substantial cosmic volume outside the 
main halo. 


N-body and hydrodynamical methods. The Ramses code 
is based on the adaptive mesh refinement technique, which 
solves the hydrodynamical equations on progressively refined 
grids. The hydrodynamical solver is based on a second-order 
Godunov scheme for ideal gases with an equation of state 
7 = 5/3. High-resolution dark matter particles have a mass 
of 1 O®M 0 , and the highest spatial resolution is physical 5 
kpc (maximal refinement level 18), with a mass-based quasi- 
Lagrangian refinement strategy. Due to the added expense of 
modelling the baryons, this first set of Rhapsody-G simula¬ 
tions has eight times lower mass resolution than the original 
Rhapsody A-body sample. 


Cooling and star formation. Our simulations follow the 
subgrid cooling model from Sutherland & Dopita (19931, im¬ 
plemented by Teyssier et al. | |2011[ | for Ramses. The star 
formation rate follows p* = e*(pgas/tff), with the star forma¬ 
tion efficiency e*=0.03 and the free-fall time tg = ■y/37r/(32Gp), 
where p is the total mass density. 

AGN feedback. We modify the AGN feedback model in 


Martizzi et al. (2012 20141, which was based on Booth & 


Schaye (20091 and Teyssier et al. (20111. In this implementation. 


supermassive black holes are modelled as sink particles, which 
grow based on mergers and Bondi-Hoyle accretion with a boost 
factor a, with an upper limit set by the Eddington rate. The 
thermal energy associated with the accretion is not released 
until the temperature reaches a certain threshold Tmin- We 
choose Q = (n_f//n*)^ when uh > n, = 0.1 P/cm? (uh is the 
gas density) and a = 1 otherwise; Tmin = 10^ K. The AGN 
thermal energy is injected into a region of four times the cell 
size. We do not implement kinetic AGN feedback associated 
with jets; it has been shown in Dubois et al. (2012 I that the 
kinetic feedback does not significantly affect the bulk gas and 
stellar mass. 


In [Teyssier et al.| ( |2011| ) and [Martizzi et al.j ( [2012| [2014[ ) , 
the feedback energy was distributed based on a volume- 
weighted approach, whereas for Rhapsody-G we adopt a mass- 
weighted approach. This implementation results in an effect 
similar to the quasar mode feedback provided by the radiation 
pressure of AGN (e.g. Debuhr et al. 20121. In addition, we 
require that black holes form at the centre of gas clumps with 
accretion rate > 3 OM 0 /yr, and we use the gas clump finding 
algorithm developed by Bleuler et al. (20151. As we will show 
in Hahn et al. (in preparation), this implementation results in 
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Main halo info (z = 0) 

ID C200 -^5OOc[M0] 

No. of well-resolved haloes 

z = 0 z = 0.5 z = 1 

572 

7.03 

5.66x1014 

2 

2 

1 

337 

5.19 

6.61x1014 

1 

2 

1 

377 

4.79 

4.91x1014 

2 

3 

3 

348 

4.62 

6.22x1014 

4 

4 

3 

653 

4.47 

3.77x1014 

3 

4 

2 

361 

4.41 

5.98x1014 

5 

5 

5 

448 

4.40 

5.81x1014 

4 

3 

1 

545 

4.40 

5.13x1014 

1 

2 

1 

211 

3.65 

5.13x1014 

5 

6 

1 

474 

3.55 

1.32x101® 

3 

6 

7 

Totals 

30 

37 

25 


Table 1. Numbers of distinct haloes in the high-resolution regions 
with Msooc > 5 X IO^^Mq. We sort the list by the concentration of 
the main halo, C200 = ^200/^s, at 2 = 0. 


Figure 1. Mass and redshift of the haloes used in this work. Each 
colour represents one zoom-in simulation centred on a main halo 
with Msooc ~ 6 X 10^‘^Mq at 2=0. We include three snapshots for 
each simulation: 2 = 0, 0.5, and 1. 


a halo mass-stellar mass relation in agreement with |Kravtsov| 


et al. (20141. 


Post processing. We modify the phase-space halo finder 
Rockstar ( [Behroozi et al.|2013[ ) to include multi-resolution, 
multi-species particles and gas. We treat the leaf-cells of the 
adaptive mesh refinement tree as pseudo-particles, thus al¬ 
lowing for direct integration of all matter components in our 
haloes. In this work, we use the dark matter density peak as 
the centre of the halo, which closely coincides with the stellar 
mass density peak in most cases. As we focus on the bulk 
cluster properties for radius > i? 2500 c, the choice of centre does 
not affect any of the results presented here. 

We adopt the same flat A cold dark matter cosmology as 
in the Rhapsody simulation. The cosmological parameters in 
the simulations are as follows: SIm = 0.25; Qa = 0.75; flf, = 
0.045; h = 0.7. Our cosmic baryon fraction value is Q.b/^M = 
0.18, which is slightly higher than the value recently reported 
by Planck (0.155, see Planck Collaboration et al.||2014|. 


2.2 Sample selection 

For this study, we include all haloes in the high-resolution 
region having Msooc > 5 x 10^^ Mq at output redshifts of 2 
= 0, 0.5, and 1. These haloes all satisfy the condition that 
the mass fraction contributed by low-resolution dark matter 
particles is below 10“^. Since each simulation encompasses a 
high-resolution sphere of 8h~^Mpc centred on the main halo 
at 2 = 0 and progressively larger high-resolution regions at 
earlier times, we are able to include 92 haloes in total. 

Fig.[T] illustrates the masses of the haloes at the three 
redshifts. Each colour represents a zoom-in simulation. For each 
simulation, the main halo and its most massive progenitors are 
shown with larger symbols, while other progenitors and nearby 
high-resolution haloes above the mass threshold are shown 
with smaller symbols. The three symbol types correspond to 
the three redshifts studied here. The symbol styles will be 
repeated in the figures below. Table lists the numbers of 


Symbol 

Quantity 

Mean (Rsooc) 

f* 

Stellar mass fraction 

0.023 

/h 

Hot gas mass fraction (kT >0.1 keV) 

0.149 

fc 

Cold gas mass fraction (kT ^0.1 keV) 

0.005 

k 

Gas mass fraction, /g = /c -h /h 

0.154 

fc* 

Galactic mass fraction, fc* = f* + fc 

0.028 

/b 

Baryon mass fraction, /b = /* + /h + fc 

0.177 


Table 2. Notation used for baryon mass components. The first three 
quantities are derived from the Ramses output; the rest are derived 
from these quantities. All fractions are relative to the total mass 
within some chosen density contrast. 


haloes derived from the outputs of each simulation, with the 
first column giving the original Rhapsody halo ID. 


2.3 Statistical error estimates 

Our main sample consists of all haloes found in the high- 
resolution regions of the 10 simulations at 2 = 0, 0.5, and 1. 
While the time spacing between these redshifts corresponds 
to several dynamical times, we do not assume statistical in¬ 
dependence among different redshifts for a given halo. When 
estimating statistical errors in quantities presented below, we 
treat the halo ensemble extracted from each simulation as 
independent. To estimate uncertainties, we jackknife resample 
using 10 degrees of freedom, eliminating one ensemble of haloes 
at each round. 


3 GALAXY CLUSTERS AS NEARLY CLOSED 
BOXES FOR BARYONS 

In the spherical collapse model of dark matter haloes, the 
tnrnaround radius sets the scale within which the cosmic mix 
of baryonic and cold dark matter should be conserved | |Gunn| 
fc Gott|1972 1. Assuming the influence of gravity and collisional 


shocks, Bertschinger (19851 developed self-similar solutions 


for both collisionless and collisional, ideal fluids, finding that 
both flnids approached similar radial profiles when expressed 
in units of the turnaround radius. Furthermore, the solution 
for the mixed case of a collisionless fluid with a minority ideal 
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f»( <I^500c) 4s( <^500c) 

Figure 2. Anti-correlation between gas and stellar mass fractions. Left: the gas mass fraction using both hot and cold phases, /g, and stellar 
mass fraction, /*, within i? 500 c exhibit an anti-correlation, r = —0.72 it 0.02. Right: the hot (fcT > 0.1 keV) gas fraction, /h, and galactic 
(cold gas and stellar) mass fraction, /c*, exhibit a stronger anti-correlation, r = —0.79 it 0.02. The symbol styles are the same as used in Fig.j^ 
The dashed line in each panel shows the perfect anti-correlation case for which the sum of the two components plotted equals the assumed 
cosmic mean baryon fraction, 0.18. 


gas component showed no radial separation; the local interior 
baryon fraction reflects the cosmic mean value at all radii. 

Taking this model to its logical extreme, let us imagine 
spherical collapse around a local perturbation consisting of 
radial shells made of cold dark matter, galactic stars, and 
smooth intergalactic gas that can shock but is unable to cool. 
Collapse of these shells would create a cluster in which the 
local interior baryon fraction was unbiased at all radii, and in 
which the mix of stars and gas interior to a given radius would 
reflect the initial values imposed on the shell layers. 

In this section, we show that this naive expectation is 
respected to a surprising degree for the case in which baryons 
experience complex astrophysical processes associated with 
galaxy formation in a fully three-dimensional, hierarchical 
clustering environment. 


3.1 Correlations among baryon mass fractions 

For each halo in the sample listed in Table we identify all 
material within a sphere that encompasses a total mass density 
contrast Ac with respect to the critical density. We measure 
the mass in cold dark matter and in all baryonic components to 
derive the total mass. We examine the six baryon mass fractions 
listed in Table Hot and cold gas phases are defined using a 
temperature cut of kT = 0.1 keV, which is approximately the 
threshold of X-ray emission. At 2 = 0, the cold gas fraction /c 
is generally very small. 

Fig. [2] shows the behaviour of different baryon mass frac¬ 
tions measured within i? 5 ooc. The point colouring and styling 
is the same as used in Fig. In each panel, the dashed line 
gives the simple expectation in which the sum of the two com¬ 
ponents plotted equals the cosmic mean, (component 

correlation coefficient of —1). 

The left panel plots gas mass fraction, /g, against stellar 
mass fraction, /*. The two components have a rank correlation 


—0.72 ± 0.02, where the error bar is calculated by jackknife 
resampling by removing one of the simulation sets in turn. In 
the right panel, we shift the mass in cold gas within i? 5 ooc 
to the stellar component, and plot the hot gas fraction, /h, 
against the total galactic fraction, /c*. This split, which more 
closely represents material inside and outside of galaxies, leads 
to a stronger anti-correlation of —0.79 ± 0.02. 

These anti-correlations are not entirely surprising, given 
that the cold gas and stellar mass originated from the cooled 
hot-phase gas. Nevertheless, we note that /g and /, differ 
by almost an order of magnitude, and that the dynamics 
within f? 5 ooc are strongly non-linear and different for collisional 
and collisionless components. Therefore, such a high degree of 
correlation is a non-trivial finding. 

Some of this anti-correlation is driven by trends in mean 
baryonic content with halo mass. Massive clusters are observed 
to have higher /g and lower /* than lower-mass systems, re¬ 
flecting a lower star formation efficiency in higher mass haloes 


(e.g. Gonzalez et al.||2007 

Giodini et al.|2009 Andreon|2010 

Zhang et al.||20111 |Lin et al.||20121 |Gonzalez et al.|20131 |Chiu| 

et al.|2014 

1 . As we will show below for our simulations, models 


with AGN feedback are capable of reproducing these trends 
(also see, e.g. McCarthy et al.|2011 Planelles et al.|2013 Le 


|Brun et al. |2C)14[). 


We fit the mean dependence of /g and /* with halo mass 
and remove these trends to examine correlations between the 
residuals, 6fg = /g — {fg\M) and Sft = f* — (/*|M). The mean 
behaviour is derived from a logarithmic fit. The correlation 
coefficients within i? 5 ooc decline slightly, to —0.63 ± 0.02 and 
—0.69 ± 0.01 for the left and right panels of Fig. respectively. 
In Table 1^ we list rank correlation coefficients at various scales 
with the mass trends removed. 


It is interesting to ask whether /g is also correlated with 
other galaxy properties. We define the stellar mass associated 
with individual galaxies as the mass within the isophotal con- 
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tour of 25 mag/arcsec^ in the V band, measured along the 
direction of the angular momentum. In our simulations, the 
total stellar mass inside Rsooc is strongly correlated with the 
stellar mass of the brightest central galaxy (BCG). Therefore, 
unsurprisingly, /g and the stellar mass of the BCG are also 
significantly anti-correlated, with a slightly weaker coefficient 
of-0.67 ±0.03. 

On the other hand, the number of galaxies, or specifically 
the ratio between galaxy number and cluster mass (A^gai/Aftot), 
has a much weaker anti-correlation with /g (r = —0.42 ± 
0.04). Here A^gai is the number of galaxies with stellar mass 
above 10 ^^Mq within Rsooc’, we have tested galaxy stellar 
mass thresholds between 10^^ and lO^^M©, and the correlation 
stays approximately constant. The main reason is that the 
lower-mass haloes In our sample are more strongly dominated 
by the BGG; that is, although they have higher stellar mass 
fraction, their A^gai/Mtot is still low. Nevertheless, if we focus 
on a narrow mass range, Ngai/Aftot is expected to correlate 
with /* and thus anti-correlate with /g. When we consider 
the most-massive haloes in each snapshot, we find a slightly 
stronger correlation (r = —0.49 ± 0.07). 

We note that the stellar mass unassociated with individ¬ 
ual galaxies, the intracluster light (ICL), is difficult to ob¬ 
serve, because it requires observations with high sensitivity 
at very low surface brightness. Here we consider the stellar 
mass associated only with the BCG and satellite galaxies, 
within the 25 mag/arcsec^ isophote mentioned above. When 
we sum all the stellar mass associated with galaxies with stel¬ 
lar mass > IO^^'^Mq inside Rsooc, we still find a significant 
anti-correlation between /g and this ICL-excluded /* (—0.69). 
This anti-correlation also weakly depends on the galaxy stellar 
mass threshold. However, we find that this anti-correlation is 
largely driven by the BCG; when we exclude the stellar mass 
of BCG, this anti-correlation is largely diminished (—0.3). 

We caution that the hydrodynamic and gravitational reso¬ 
lution of ~ 5 kpc in our simulations is insufhcient to reliably 
model the ICL component. Using the 25 mag/arcsec^ threshold, 
we find an average ICL mass fraction of 55% of the stellar mass. 
This fraction is generally higher than observed, reflecting our 
models’ inability to resolve the half-light radii of all but the 
brightest galaxies. A simulation suitable for studying ICL will 
require higher spatial resolution. Nevertheless, the consistency 
of the gas-stellar correlation coefficients derived with and with¬ 
out the ICL indicates that our results are not being driven by 
the ICL component. 


3.2 Scale dependence 

We expand the results shown in Fig. [^to explore the scale 
dependence of baryon component correlations. We compute 
statistics within radii that encompass overdensities of Ac = 
2500, 500, 200, 50, and 10. For the two lowest overdensities, we 
use only the most massive progenitor halo in each simulation 
to avoid contamination from low-resolution particles. 

Fig. [3] shows the correlation coefficients as a function 
of density contrast (large radii are towards the right). The 
circular symbols give the correlation between /g and /*, while 
the triangular symbols correspond to the correlation between 
/h and /„. The left panel gives the raw correlation before 
removing the mass trends, while the right panel removes the 
effect of the mass trends. In Table we list values of the latter. 

The anti-correlation between hot phase and galactic (cold 
gas and stars) components is stronger than that between gas 


and stars at all radii. In both panels, the former Is very close 
to —1 at Ac = 10. At low density contrasts, the Influence of 
trends with halo mass is weak. At higher density contrasts, 
the trends with halo mass are more important, and correcting 
for them dilutes the raw correlations by ~ 0.1. Still, within 
772500 c the correlation between hot and total galactic baryon 
residuals is —0.71 ± 0.02. This strong covariance indicates that 
the simple closed-box model remains approximately valid even 
at radii where complex galaxy formation physics is at play. 


3.3 Dependence on the astrophysical treatment 


We caution that the anti-correlation between /g and /« pre¬ 
sented above is based on one particular astrophysical treatment 
of star formation and feedback. To test whether these results 
are sensitive to this treatment, we employ a set of 51 cluster 
simulations from |Martizzi et al.| ( |2014| M14 hereafter). These 
simulations are based on the same code and methods as used 
here, but they employ a volume-weighted AGN energy injection 
model as opposed to the mass-weighted model in Rhapsody-G 
(see Section]^. With the volume-weighted feedback, energy 
injection within the core is more efficient and the BCG star for¬ 
mation is more strongly suppressed compared to mass-weighted 
feedback. 

The M14 sample reproduces many of the trends presented 
in this paper, but there are some differences. In particular, 
the stellar mass fraction, while similar in the mean to that 
of Rhapsody- G, has a factor of 2 smaller scatter at a fixed 
halo mass (17% in M14 compared to 34% in Rhapsody-G at 
Rsooc)- The observed value reported by [Kravtsov et ah (20141 
of 0.29 ± 0.09 (fractional scatter, not dex) slightly favours the 
latter but is marginally consistent with both. 

Along with the smaller scatter in /* at R^ooc the anti¬ 
correlation between /g and /* at a given mass in the M14 
sample Is also reduced by a factor of 2, to r = —0.35. The 
anti-correlation thus persists qualitatively, but the M14 sample 
behaves somewhat less like ‘closed boxes’ within their non¬ 
linear regions than do the Rhapsody-G sample. 

The small /* scatter in the M14 sample is associated with 
the implementation of the AGN feedback. The exact value of 
the correlation coefficient between /g and /* thus depends on 
how AGN feedback is modelled. Uncertainties in predicting 
stellar masses in simulations (e.g. Ragone-Figueroa et al.|2013 


Martizzi et al. 20141 imply concurrent uncertainties in the 
correlation between /g and /*. Measuring the anti-correlation 
between hot gas and galactic mass fractions, discussed in Sec¬ 
tion]^ thus offers a means to constrain details of the feedback 
model. 


4 BARYON MASS AS A LOW-SCATTER 
PROXY OF TOTAL MASS 

The anti-correlations in baryonic components shown above 
reflect the fact that the sum of these components has lower 
scatter with respect to total halo mass than does each compo¬ 
nent. Fig. shows scaling relations for gas mass, stellar mass, 
and baryon mass fractions as a function of total halo mass, 
measured within Rsooc- The left panel shows that fg is very 
tightly correlated with the total mass, with a scatter of 8% 
and a slope of 0.08. We note that this scatter is slightly lower 
than previously reported using a wider range of halo masses; 
e.g. [Kravtsov et ah] (|2006| reported a fractional scatter of 0.107 
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Figure 3. Anti-correlations between gas and stellar mass fractions as a function of overdensity Ac- Small Ac (large radius) is to the right. At 
large radius (low Ac), the anti-correlation is significantly stronger, but it remains robust at small radius. The anti-correlation is stronger 
between Sf^iSfc* (green) than between (5/g—<5/* (blue); 5f is defined as / — (/|M), i.e., with the mass dependence subtracted. 





^500c[^o] ^SOOcl^ol ^500c[^o] 

Figure 4. Scaling relations of gas mass fraction (left), stellar mass fraction (middle), and their sum (right) with total halo mass at 2 = 0, 0.5, 
and 1. The symbol styles are the same as used in Fig.^ The general trends with mass are consistent with those observed, while the exact 
values of the slopes may differ from those observed because our sample is incomplete at the low mass end. 


Summary of correlation and scatter for gas, stellar, and baryon mass fractions 



Rank correlation 


Fractional scatter 


(■5/g. <S/*) 

(<5/h, SM 

/g 

f* 

/b 

R 2500 C 

-0.63 ±0.02 

-0.69 ±0.01 

0.19±0.005 

0.38±0.007 

0.11±0.002 

R 500 C 

-0.66 ±0.02 

-0.69 ±0.02 

0.08±0.002 

0.34±0.005 

0.047±0.002 

R 200 C 

-0.66 ±0.02 

-0.75 ± 0.02 

0.062±0.002 

0.32±0.005 

0.038±0.001 

R 50 C 

-0.92 ± 0.02 

-0.92 ± 0.01 

0.024±0.002 

0.26±0.02 

0.01±0.0007 

Rioc 

-0.89 ±0.03 

-0.980 ± 0.003 

0.015±0.001 

0.22±0.02 

0.0047±0.0003 


Table 3. Correlation coefficients between gas and stellar mass fraction (second column), between hot gas and galactic baryon mass fraction 
(third column), as well as the fractional scatter in gas mass (fourth column), stellar mass (fifth column), and baryon mass (sixth column) 
fractions, at radii characterized by different overdensities. 
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Figure 5. Scatter in baryon mass proxies as a function of radius 
defined by the enclosed overdensity, A^. The relation has 

the lowest scatter at all radii, falling to 0.5% at Ac = 10. 


(based on R^ooc, see their table 2). On the other hand, our 
scatter is similar to the observed value recently reported by 
|Mantz et al.| ( |2014[ ) . For a sample of 40 relaxed clusters ob¬ 
served with Chandra, they have found 7% scatter in /g within 
a shell of 0.8-1.2i?25ooc. 

The middle panel shows that /* has a much larger scat¬ 
ter of 34% (0.13 dex) and a slope of —0.21. In observations, 
Kravtsov et al. ( 2014|) recently reported that the Mt-Msooc 


relation has a similar scatter (0.11 ±0.03 dex). The right panel 
shows that when /g and /* are combined to form /b, the scat¬ 
ter is 4.7%, which is smaller than using either /g or /* alone. 
The overall trends with mass reflect that low-mass haloes have 
higher efficiency of turning gas into stars, but at the same time 
they tend to have larger baryon depletion. 

In Fig. we show the fractional scatter of the different 
baryon mass components as a function of scale. The fractional 
scatter of M, is the largest. The scatter In Mg is smaller and 
declines more rapidly with increasing radius. The fractional 
scatter in the overall baryon mass, Mb, is typically a factor of 2 
smaller than that of the gas mass, declining to 0.5% at Ac = 10. 
Overall, the baryon fraction scatter scales approximately as 
^ 0.018(Ac/100)° ®». 

The reduced scatter in the total baryon mass can also be 
understood in the context of joint property selection. Under the 
assumption of a joint lognormal distribution for two observables 
at a fixed halo mass, the scatter of halo mass obtained by joint 
selection in these observables is given by ( Evrard et al.|20i4 1 

_2 1-R 

—2 , —2 r , -1 —1 ’ 


-2 I -2 o -1 -1 ’ 

<^1 + <^2 “ 2ra, On 


where ui and 02 are the fractional scatter in the two observables 
at a fixed mass, r is their correlation coefficient, and (Tjoint is the 
resultant scatter in halo mass under joint selection. Evaluating 
this expression using cri = 0.08, (J 2 = 0.34, and r = —0.72, we 
obtain ujoint = 0.047, which matches the fractional scatter we 
have found in Mb. 


5 OBSERVATIONAL IMPLICATIONS 

In this work, we predict a significant anti-correlation between 
gas and stellar mass fractions in the virial regions of galaxy 
clusters. Here, we discuss a few possible approaches that could 
be applied to current and future observations to test for this 
signature. The fundamental challenge lies in understanding 
sources of noise that would drive the measured covariance away 
from the intrinsic value. 


5.1 Individual measurements of /g and /* 


A direct approach is to measure gas, stellar, and total masses 
for an ensemble of clusters and examine the correlations in 
the component fractions. Ideally, the noise in such measure¬ 
ments would be smaller than the intrinsic scatter, which is 
approximately 10% for Mg and 30% for M,. 

Gas masses can be determined fairly reliably; |Rozo et al.| 
( |2014[ | showed that, after accounting for aperture biases driven 
by total mass errors, gas mass estimates of individual clus¬ 
ters deviate in the mean by just a few percent across three 
independent observing teams using different X-ray telescopes 
and analysis methods. Statistical uncertainties from photon 
statistics can be made small for estimates at Ac ^ 500 in low 
redshift clusters. 

Deriving the stellar mass of clusters from optical images 
is more difficult. Stellar masses of bright galaxies are sensitive 
to the surface brightness profiles ht to the photometry (e.g. 
Bernard! et al.|[2013l, and the revised fits by Kravtsov et al. 


(2014) pushed BCG stellar masses of SDSS galaxies up by 


factors of 2-4. Systematic uncertainties in cluster membership 
assignment can produce errors in the stellar mass contributed 
by satellite galaxies (e.g. Lin et aL]|2004 1 , and uncertainties in 
the stellar population synthesis models can also contribute (e.g. 


Conroy et al. 20091. Finally, measuring and defining intracluster 
light involve additional uncertainties (e.g. Lin fc Mohr|200'4 


[Gonzalez et al.|2007[ [Budzynski et al.|2014|). 

These additional uncertainties would dilute the measured 
correlation. Adding random fractional error, ao, to the stellar 
mass estimates alone would decrease the correlation coefficient 
by a factor of ct* / ^/o^g + where cr* is the intrinsic scatter of 
roughly 30% (Table |^. For (tq = u*, the correlation coefficient 
would reduce to ~ —0.5. 

Statistical error of a few tens of percent is expected in 
total mass estimates of individual clusters derived from either 


et al.|2007 Rasia et al.|2012 2014 Nelson et al.|2014 

) or weak 

gravitational lensing (e.g. 

Becker & Kravtsov||2011 

Oguri & 

Hamana||2011| 

Bahe et al. 

|2012||Rasia et al.||2012||. However, 

Donahue et al. 

(2014) have provided evidence that combined 


strong and weak gravitational lensing models reduce the scatter 
in total mass estimates. Since error in total mass induces a shift 
of the same sign for /g and /*, such errors would smear out 
the trend seen in Fig. by introducing scatter perpendicular 
to the anti-correlation. A potential alternative is to avoid total 
mass estimates and select instead on a low-scatter mass proxy, 
such as X-ray temperature or the product between gas mass 
and X-ray temperature (lx). One could then examine how 
Mg and M* within fixed metric apertures covary within that 
sample. This approach would take advantage of the relatively 
weak radial dependence of the correlations exhibited in Fig. 

Fig.§ compares our results with several observed clusters 
in the literature. Our simulations are represented by the grey 
histogram in the background (the same data points as in 
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Figure 6. Comparison of component mass fractions in our simu¬ 
lations (grey 2D histogram) with observed clusters from [Lagana] 
|et al.| ( (2011[ | and |Gonzalez et al.| | |^13| l. The component fractions 
are normalized to the universal baryon fraction, fb^univ = 
which for our simulations is 0.18. For the observations, we use the 
recent Planck satellite value of 0.155. 


the left panel of Fig. [^. We include observed clusters with 
Msooc > 5 X IO^^Mq published in Lagana et al. (2011) and 
[Gonzalez et al.| ( |2013| |. We note that these data sets are based 
on different mass calibration techniques. Lagana et al. (2011) 
used a Mg-Msooc relation to calculate Msooc and used the 
Schechter function for [Gonzalez et al.| ( |201^ used a Tx- 
Msooc relation to calculate M^ooc and carefully accounted for 
ICL when calculating M*. 

In our simulations, Qb/^m = 0.18, which is higher than 
the value of 0.155 recently constrained by the Planck satellite 
( Planck Collaboration et al.[[2014 ). To account for this differ¬ 
ence, we normalize the simulated component fractions by our 
flb/^m, and the observed fractions by the Planck value. As 
shown in Fig.[^ the statistical error bars on /g and /* are fairly 
large, and there are systematic effects that we cannot correct 
for. Setting these caveats aside, the simulated and observed 
trends in Fig. are roughly consistent. 


5.2 Multi-property scaling relations 


Since directly measuring the correlation between /g and /* 
is challenging, a practical alternative is to consider statistical 
effects on scaling relations for property-selected samples. The 
basic idea is as follows: if we select a subset of clusters based on 
property Sa, then the mean and scatter of a second property, 
Sb, will depend on how Sa and Sb are correlated at a given 
mass. Following the formalism of Evrard et al. (2014), we 
assume a power-law mass-observable scaling, so that 

{Sa\p) = T^a + OafJ., ( 2 ) 


where Sa = InS'a, fJ, = InM, and tto gives the normalization in 
the chosen units. Denoting the scatter in Sa at a given mass 
by (Ta, then the scatter in /r at a given Sa is cr^|a = Oalo-a- 


The mass scaling of the second observable, Sb, follows 
similar notation. Let Tab be the correlation coefficient of Sa 
and Sb at a fixed mass. Finally, the convolution from mass to 
observed signal requires knowledge of the mass function, which 
can be approximated as n{M) = Ae~^'^, with —/3 the local 
logarithmic slope of the mass function. 

For a sample selected on observable Sa, the mean value of 
Sf, is given by 


-)- Oib 


{p\Sa)+Prab 


(3) 


where (/r| Sa) is the mean halo mass selected by Sa - The variance 
of Sb is given by 

^b\a ^6 [^/i|a T 2 Tab (4) 


These expressions show that if Tab < 0, then the mean of Sb 
will be biased low, and its scatter will be larger compared to 
the case of no correlation. 

The effect on the mean will generally be small. For example, 
if the two observables have 20% mass scatter, and /3 ~ 2.5, then 
the hnal term in Equationis a shift of O-lCab in InM. Such 
small shifts, below 10 percent in mass, are currently challenging 
to measure, since the systematic errors in mass are of similar 
or larger magnitude. In addition, measuring this shift requires 
accurate knowledge of the observable-mass normalization, TVb, 
as well as the mean selected mass. 

Analysing the variance is a potentially simpler alternative. 
If properties a and b have comparable mass scatter, then all 
the terms on the right-hand side of Equation will be of the 
same order. As an example, we consider recent observational 
results involving galaxy richness (as Sa) and gas mass {Sb) in 


Rozo & Rykoff (2014), which are summarized in their table 2. 


The mass scatter at a fixed richness is approximately cr^|Q = 
0.25 ( [Rykoff et al.|[^012[ ), while the mass scatter at a fixed 
Mg is approximately a^\b = 0.1 (Mantz et al. 2010| 2014). 
Based on matching existing X-ray data to the optically selected 
redMaPPer cluster sample, [Rozo fc Rykoff[ ( [2014[ ) report a 
scatter in gas mass at a fixed galaxy richness of (Ti,|a = 0.212 ± 
0.032, and report a slope for the M^-Mx relation of ab = 
0.72 ±0.12. 


Evaluating Equation|^with these values gives Tab = —0.28, 
a slight hint of an anti-correlation. However, the exact value 
of Cab sensitively depends on the slope ab, which is poorly 
constrained. Decreasing the slope by its one sigma uncertainty, 
to ab = 0.65, leads to an estimate of Tab = —0.68. Similar 
exercises based on larger samples of homogeneously determined 
mass estimates will offer a more robust means to test for non¬ 
zero covariance in stellar and hot gas mass fractions. 

A related test is to probe the variance in total halo mass 
under joint property selection (see Equation]^. Using optical 
and X-ray samples for which both M, and Mg are accurately 
measured, one could first use lensing total masses to estimate 
the scatter in Ma for each property selection. Fitting the 
lensing masses in the joint selection of both properties would 
provide a fundamental plane with variance reduced by an 
amount given by Equation Selection effects would need to 
be carefully modelled in such a study. 


6 SUMMARY 

We present an analysis of the various baryonic mass components 
- stars, hot and cold gas - in a sample of ~ 100 massive haloes 
derived from the Rhapsody-G cosmological hydrodynamic 
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simulations. These simulations include state-of-the-art models 
for gas cooling and star formation, as well as for energy injection 
through supernovae and AGN. Our findings can be summarized 
as follows: 


• At a fixed total halo mass, stellar and gas mass fractions 
are significantly anti-correlated in the non-linear regions of 
haloes, with r = —0.66 ± 0.02 at flsooc This correlation is 
further enhanced if we split gas into hot and cold gas, and 
correlate the hot gas fraction with the cold gas plus stellar 
mass fraction (r = —0.69 at /?5ooc). 

• Due to this anti-correlation, total baryon mass has a scat¬ 
ter with respect to the total halo mass that is lower than either 
gas mass or stellar mass. At 7?5ooc, the baryon mass has ap¬ 
proximately 5% scatter, suggesting that joint cluster selection 
using accurate gas and stellar mass measurements can achieve 
up to 5% selection in total mass. 

• With increasing radius, the anti-correlation between /g 
and /* approaches -1, the closed box expectation, and the 
baryon mass scatter declines to 0.5% at Ac = 10. 


It is currently challenging to accurately measure the anti¬ 
correlation between /g and /* in observations. Scaling laws 
with well measured slopes, intercepts, and standard deviations 
are required for large samples. To obtain empirical constraints 
on this correlation in massive clusters, joint survey studies 
that combine gas mass and stellar mass selection with lensing 
masses and/or additional independent mass proxies (from X-ray 
temperatures, Yx, caustic masses, galaxy velocity dispersions, 
etc.) are needed. 

Comparison to Ramses simulations that use a different 
AGN feedback scheme indicates that the results are qualita¬ 
tively robust but quantitatively dependent on the feedback 
method. In light of the different simulation results from differ¬ 


ent implementations of AGN feedback (e.g. Ragone-Figueroa 
et al.|[MT3 Martizzi et al.|[2014l and between adaptive mesh 


refinement and smoothed-particle hydrodynamics methods (e.g. 
|Frenk et al.f|1999j |Rasia et al.||2014{ [Sembolini et al.||2015P , it 

is important to address the anti-correlation found here using 
different simulation techniques, more detailed physical models, 
and different subgrid models for feedback processes. The closed- 
box result must hold at sufficiently large radii, but the detailed 
scale dependence of the covariance in baryon components is 
likely to exhibit model-dependent features. 
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